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1. INTRODUCTION 

Image processing applications plays a vital role in the present digital world. In that satellite imaging 
has many useful applications in agriculture, remote sensing, urban planning and mineral exploration. Despite 
many useful applications, it is critical to transmit the satellite image data into the ground station in the limited 
transmission bandwidth. The satellite will be in contact with the ground station only at certain intervals. The 
data need to be stored on-board until the next communication window of the satellite. For all the reasons, it is 
very important to compress the satellite image data to make the communication and storage feasible. Image 
compression algorithms are broadly classified into lossless and lossy image compression algorithms. The 
satellite image is not only used for viewing purposes but also used to make important exploratory analysis. 
Therefore, lossless image compression [1] is the modality used in most of the applications. But in certain 
cases, to improve the image acquisition rates and to improve the compression rate lossy satellite image 
compression algorithms are made use in missions like SPOT 4 [2]. 

Nyquist sampling theorem [3] is the well-known method for signal acquisition and faithful recovery 
of signal from the measurements. Advances in compressive sensing-based methodologies gives satisfactory 
proof on retrieving the signal from a smaller number of samples when the signal is compressible in nature. 
The modality of compressive sensing is very useful when it is not practical to acquire many samples and 
when it is very expensive to collect large number of samples. In both cases compressive sensing can 
contribute by making the image acquisition hardware simpler, reducing the image acquisition time, and 
thereby achieving image compression. 
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Conventional image compression algorithms rely on the Nyquist sampling theorem where large 
numbers of samples are collected only to be discarded at the stage of compression. Conventional image 
compression algorithms in the literature are generally based on transform coding techniques. Where, the 
images are transformed using discrete cosine transform (DCT) or discrete wavelet transform (DWT). The 
coefficients are then encoded using encoding techniques like run-length encoding, Golomb coding [4], [5], 
Arithmetic coding [6], and rice coding [7]. Many conventional image compression algorithms have been 
proposed in the literature for the compression of satellite images. The consultative committee for space data 
systems (CCSDS) document [8] describes the lossless compression methods for satellite images. The image 
compression algorithm based on DWT called set partition in hierarchical trees (SPIHT), where the significant 
DWT coefficients are encoded first for the fast and efficient coding of images is presented in [9]. Aranki et 
al. [10], the field programmable gate arrays (FPGA) implementation of the lifting scheme of DWT is 
proposed. The research work in [10] and [11] discuss the concept of finding the DWT coefficients of an 
image segment, in a pipeline fashion, instead of waiting for the entire image from the push-broom scanning, 
since the values of DWT coefficients depends on the pixels in the neighbourhood. JPEG-LS algorithm for the 
lossless compression of images is presented in [12]. A low complexity context based lossless image 
compression algorithm is presented in [13]. The JPEG-2000 algorithm, which is a DWT version of DCT 
based JPEG to achieve different degrees of image quality based on different applications is proposed by the 
joint photographic experts’ group in [14]-[16]. Compressive sensing [17]-[19] is based on the principle that if 
the signal is compressible in nature, it is possible to reconstruct the signal from very few measurements. The 
compressive samples are generated using the random sampling matrix which captures the signal in very few 
measurements. 

Figure 1 shows the general block diagram of the conventional image compression algorithms like 
JPEG and CCSDS standard. As it is evident from the block diagram, the satellite image goes through many 
processing steps at the encoder stage. First the image is transformed using one of the transform coding 
techniques like DCT or DWT. The insignificant coefficients are discarded after the transform coding. The 
remaining coefficients are subjected to quantization in the next stage. The quantized coefficients are further 
encoded using entropy coding techniques like Huffman coding, and arithmetic coding. The encoder is very 
simple in case of CS based compression algorithms as shown in Figure 2. The only processing step at the 
encoder is generation of compressive sensing samples from the random compressive sensing matrix. The 
decrease in computational steps in the encoder is of great significance in case of remote sensing application 
because memory and power are the two important critical resources in case of any remote sensing 
application. 
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Figure 1. General block diagram of the conventional satellite image compression algorithm 


The most significant aspect of using CS for the compression in remote sensing applications can be 
attributed to the property of making the encoder very simple by removing all the processing steps in the 
encoder stage of the conventional compression algorithm. The compressive sensing samples are 
communicated to the ground station in case of satellite imaging applications. The disadvantage of CS based 
algorithms for compression is that the reconstruction process from compressive sensing samples is quite 
computationally expensive. This aspect can be overlooked as the reconstruction process happens in the 
ground station where resources are available in abundance. 
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The research work makes use of compressive sensing techniques to compress the satellite images. 
Compressive sensing samples are generated using a CS matrix containing random Fourier transform 
coefficients. The satellite image is compressed at various sampling rates to analyse the performance of the 
algorithm at different levels of compression. The satellite image is reconstructed from the compressive 
sensing samples by a class of interior point methods called log barrier method. The iterative technique 
minimizes the objective function with total variation regularization to solve the problem of image 
reconstruction from compressive sensing samples. The rest of the article is organized as follows; Section 2 
describes the compressive sensing technique proposed to reconstruct the satellite image from compressive 
sensing samples. Section 3 presents the experimental results using different datasets. The performance of the 
algorithm is analysed using quantitative metrics like peak signal to noise ratio (PSNR) and structural 
similarity index measure (SSIM). Section 4 concludes the research work presented in the article. 
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Figure 2. General block diagram of compressive sensing-based algorithm 
2. CS-CR ALGORITHM 
Consider an image I, which has dimensions n X n The total number of pixels in I is denoted by N. 
The image I is vectorized as x := Vec(I), where the length of the vector x is N. The results in papers 
[20]-[22] demonstrates the possibility of recovering signals from fewer number of measurements i.e., 
z = Cx ER”,P KN with considerable accuracy. The basis pursuit problem to solve for optimal value of x 
from fewer compressive sensing samples z is given by (1): 


min,||x||, such that Cx =z (1) 


The optimization problem solves for x with minimum l, norm, which explains the compressive sensing 
measurements Z. Another strategy to solve for x is given by (2): 


min, TV(x) such that Cx =z (2) 


where, TV (x) represents the total variation [23] of x. TV (x) is defined as (3), 


TV(x) = Les (DD)? + (D?;D” (3) 
where, 
lii> hl; i<n 
DEI = { it1 Lj 4 
a 0 i=n a 
Ti ja, li; j<n 
vi. 1,j+1 LJ 
De fb ae (5) 


The problem in (2) can be generalized as (6), 
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min, <k,,X> suchthatCx=z,F;(x)< 0 i=1...m (6) 


where, F; represents the m linearity constraints. The optimization problem in (6) can be converted into a 
series of sub problems with linearity constraints using the Log-barrier method [24] as (7): 


min, < kọ, x > EY: -log(-F(x)) such that Cx = z (7) 


Log-barrier method augments the inequality constraints into the objective by log penalty function. 
At each iteration L, the value of t is increased. When Tt is sufficiently large, the solution x" approaches the 
optimal value of the original problem in (6). The second order approximation of (7) at Ax is given by (8): 


Fo(x + Ax) = x+< Gr, Ax > +: < H,Ax,Ax > (8) 
where, Gr, represents the gradient, which is given by (9): 
1 1 
Gr, = ko + zhi ea Vi) (9) 
H,, Represents the Hessian matrix, which is given by (10): 


1 i T 1 1 
A, = 7 Linge Vie) (VFi(x)) + site A (10) 
If x is in feasible set, Ax which represents the minima of second order approximation in (8) is given by the 
following set of linear: 


T 
ale F ) (>) --16r, (11) 
p is the Lagrange multiplier for the equality constraints. 
The flowchart of the CS-CR algorithm based on Log-barrier method for the reconstruction of image 
from compressive sensing measurements is given in Figure 3. The algorithm is initialized with x°, tol, L, T, a 
in step 1. In step 2 (7) is solved using Newton method using x‘~? as initial point. In step 3, if z < tol, 
L 


algorithm is terminated, else the algorithm goes back to step 2. 


Initialize x”, tol, L, T, a 


Solve (7), using x77" as 
initial point 


Terminate 


output the value of 
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Figure 3. Flowchart of the CS-CR algorithm 
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3. RESULTS AND DISCUSSION 

Experiments are conducted using MATLB R2016b in windows 10 laptop, 64-bit operating system 
with Intel(R) 17 processor, 16 GB RAM. For the experimentation two datasets [25] have been considered. A 
band from coastal bl data of Landsat NASA is considered as the first data. The image is 1024x1024 in size. 
Another data considered for experimentation is spot la b3 from SPOT 3 CNES. The image is 500x500 in 
size. Three algorithms [26] are considered to compare and validate the performance of the CS-based 
reconstruction algorithm. The algorithm proposed in CCSDS document, uses float DWT to achieve lossy 
compression of satellite images, where the parameters are tuned to achieve different rate distortion 
performance and full frame compression. The JPEG-2000 algorithm based on float DWT in both scan mode 
and frame mode is considered for the comparison where rate allocation optimization is done using the 
Lagrange multipliers. Another algorithm considered for comparison is the set partitioning in hierarchical 
trees (SPIHT), which is based on float DWT to down sample the image into different decomposition levels. 
The methods are compared in terms of quantitative metrics like PSNR and SSIM. PSNR provides the 
information about the average numerical similarity between the original and reconstructed image. SSIM 
provides the measure of average spatial structural similarity between original and reconstructed image. The 
measures PSNR and SSIM are defined as follows: 
- Peak signal to noise ratio (PSNR): 


max(X) 


PSNR(X,X) = 20l0g:0 Man (12) 
X and X are the original image and the reconstructed image, respectively. 
- Mean square error (MSE): 
MSE(Xx,2) ==||x — 2||.” (13) 
NOON 2 
- Structural similarity index measure (SSIM): 
SSIM(X, 2) Pt: (24x03 FC1)(20x374K2) (14) 


(u%+u3+C1)(0%+03+K2) 


Hy and pg refers to the mean values of original image and reconstructed image, respectively. oy and ogrefers 
to the standard deviations of the original and reconstructed image, respectively. K, and Kzare the constants 
related to the pixel values dynamic range. Higher PSNR and SSIM is the indication for increase in 
reconstruction accuracy. PSNR and SSIM are evaluated at different sampling rates for the CS-based 
algorithm. Sampling rate is the proportion of samples extracted from the original image with respect to total 
number of samples in the original image. 

Figure 4 shows the reconstruction results using CS-based algorithm for coastal bl Landsat data at 
various sampling rates ranging from 0.1 to 0.5. Similarly Figure 5 shows the reconstruction results using CS- 
based algorithm for spot_la_b3 data at various sampling rates ranging from 0.1 to 0.5. Figure 6(a) provides 
the comparison of PSNR values for two datasets spot_la_b3 and coastal bl using conventional compression 
techniques like JPEG-2000, SPIHT and CCSDS and the proposed CS- based technique. It is evident from the 
graph that, PSNR achieved using CS-based algorithm is much higher in case of both the datasets. 

The average PSNR value using JPEG-2000 algorithm for the spot_la_b3 dataset is 42.21. The 
average PSNR value for the same dataset using CCSDS based algorithm is 41.84. Similarly, the PSNR value 
using SPIHT algorithm is given by 41.68. The PSNR value achieved using the proposed compressive 
sensing-based algorithm is 74.82 for the spot_la_b3 dataset, which is 77%, 78%, 79% better compared to 
JPEG-2000, CCSDS, and SPIHT algorithms, respectively. The average PSNR value using JPEG-2000 
algorithm for the coastal bl dataset is 51.39. The average PSNR value for the same dataset using CCSDS 
based algorithm is 51.05. Similarly, the PSNR value using SPIHT algorithm is given by 50.64. The PSNR 
value achieved using the proposed compressive sensing-based algorithm is 93.67 for the coastal bl dataset, 
which is 82.27%, 83.48%, 84.97% better compared to JPEG-2000, CCSDS, and SPIHT algorithms, 
respectively. Figure 6(b) provides the PSNR values using CS based compression technique at different 
sampling rates. It is evident that there is a steady increase in the PSNR values with increase in the sampling 
rate for both datasets. The PSNR at 0.1 sampling rate for spot_la_b3 data is 72.68. The PSNR at 0.1 sampling 
rate for coastal bl data is 90.97. Figure 6(c) provides the structural similarity index measure for the two 
datasets considered at various sampling rates using CS-based compression technique. Average SSIM value in 
case of coastal bl data is 0.9984. Average SSIM value in case of spot la b3 data is 0.9961. Figure 6(d) 
provides the run time of the CS-based compression technique at various sampling rates. The average run time 
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for the spot la b3 data is 341.97s. The average run time for the coastal bl data is 1027s. Decrease in run 
time with increase in the sampling rate indicates that the algorithm converges faster when a greater number 
of samples are available as input. Figure 6(e) and Figure 6(f) provides the compression size in bytes of the 
original and reconstructed image at various sampling rates using CS-based compression technique for the 


datasets coastal bl and spot la b3, respectively. 


Reconstruction at 0.1 sampling rate 


Reconstruction at 0.2 sampling rate 


Reconstruction at 0.3 sampling rate 


Reconstruction at 0.4 sampling rate 


Reconstruction at 0.5 sampling rate 


Figure 4. Reconstructed image at 0.1, 0.2, 0.3, 0.4, 0.5 sampling rates respectively, coastal bl Landsat data 
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Figure 5. Reconstructed image at 0.1, 0.2, 0.3, 0.4, 0.5 sampling rates respectively, spot la b3 CNES data 
[25] (All rights reserved O CNES (2007)) 
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Figure 6. These figures are; (a) PSNR using different compression techniques, (b) PSNR at various sampling 
rates using CS based compression technique, (c) SSIM at various sampling rates using CS based compression 
technique, (d) run time in seconds at various sampling rates using CS based compression technique, 

(e) original and compressed size of coastal bl image at different sampling rates, (f) original and compressed 
size of spot_la_b3 CNES image at different sampling rates 


4. CONCLUSION 

The main goal of this research work is to demonstrate the usefulness of the compressive sensing- 
based compression techniques in satellite image compression domain. The experimental results and analysis 
demonstrate the supremacy of the compressive sensing-based technique over the conventional image 
compression techniques both in terms of the reconstruction accuracy and compression rate. The compressive 
sensing techniques will be extremely useful in compression of huge hyperspectral and multispectral satellite 
image datasets as they greatly alleviate the storage and transmission bandwidth bottlenecks in remote sensing 
applications. 
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